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Abstract 

In the studies of dynamics of pathogens and their interactions with a host immune system, an 
important role is played by the structure of antigenic variants associated with a pathogen. Using 
the example of a model of antigenic variation in malaria, we show how many of the observed 
dynamical regimes can be explained in terms of the symmetry of interactions between different 
antigenic variants. The results of this analysis are quite generic, and have wider implications for 
understanding the dynamics of immune escape of other parasites, as well as for the dynamics of 
multi-strain diseases. 

1 Introduction 

In the course of evolution, pathogens have developed various methods of evading the immune system 
of their hosts. Whilst there are many contributing factors that determine individuals aspects of 
host-parasite interactions, from a more general perspectives strategies of immune escape can be 
divided into two major classes. In the first class, parasites remain largely invisible to the immune 
system of their host through an extended period latency when production of new viruses or bacteria 
inside the host organism is very small or absent. In this case, the pathogens do not trigger immune 
response, and thus are able to remain undetected in their hosts for long periods of time. Notable 
examples of such pathogens include various viruses, such as members of herpes virus family |39| and 
retroviruses [U H2] . Whilst this method of immune escape is largely unavailable to bacteria due to 
fundamental differences in replication strategies, several types of bacteria have shown persistence for 
long periods of time with little evident replication; examples include mycobacteria and T. Pallidum 
causing syphilis [50] , 

Another possible strategy of immune escape is that, in which pathogen is actively replicating in 
its host, but this replication is dynamically regulated by the host immune system. There are two 
ways how this regulation can be achieved: the pathogen can either reach a certain chronic state 
(equilibrium) through the balance of its proliferation and destruction (see |43] for Trypanosoma 
cruzi example), or it can go through the process of antigenic variation, whereby it keeps escaping 
immune response by constantly changing its surface proteins, thus going through a large number of 
antigenic variants. Perhaps, the best studied pathogens relying on this strategy for immune escape 
are Plasmodium falciparum and African Trypanosoma exemplified by Trypanosoma brucei, with 
other examples including several families of viruses, bacteria and even fungi [HI H7J I33J [36j l38j 163]. 
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In the case of T. brucei, the organism that causes sleeping sickness, parasite covers itself with 
a dense homogeneous coat of variant surface glycoprotein (VSG). Genome of T. brucei has over 
1000 genes that control the expression of VSG protein, and switching between them provides the 
mechanism of antigenic variation |41j . What makes T. brucei unique is the fact that unlike other 
pathogens, whose antigenic variation is typically mediated by DNA rearrangements or transcrip- 
tional regulation, activation of VSGs requires recombination of VSG genes into an expression site 
(ES), which consists of a single vsg gene flanked by an upstream array of 70 base pair repeats and 
expression site associated genes (ESAGs). T. brucei expresses one VSG at any given time, and the 
active VSG can either be selected by activation of a previously silent ES (and there are up to 20 ES 
sites), or by recombination of a VSG sequence into the active ES. The precise mechanism of VSG 
switching has not been completely identified yet, but it has been suggested that the ordered appear- 
ance of different VSG variants is controlled by differential activation rates and density-dependent 
parasite differentiation |41l [60] . 

For the malaria agent P. falciparum, the main target of immune response is Plasmodium falci- 
parum erythrocyte membrane protein- 1 (PfEMPl), which is expressed from a diverse family of var 
genes, and each parasite genome contains approximately 60 var genes encoding different PfEMPl 
variants |25j. The var genes are expressed sequentially in a mutually exclusive manner, and this 
switching between expression of different var gene leads to the presentation of different variant sur- 
face antigens (VSA) on the surface of infected erythrocyte, thus providing a mechanism of antigenic 
variation |11| |4~9] . In all cases of antigenic variation, host immune system has to go through a large 
repertoire of antigenic variants, and this provides parasites with enough time to get transmitted to 
another host or cause a subsequent infection with a different antigenic variant in the same host. 

Despite individual differences in the molecular implementation of antigenic variation, such as, 
gene conversion, site-specific DNA inversions, hypermutation etc., there are several features com- 
mon to the dynamics of antigenic variation in all pathogens. These include ordered and often 
sequential appearance of parasitemia peaks corresponding to different antigenic variants, as well 
as certain degree of cross-reactivity. Several mathematical models have been put forward that aim 
to explain various aspects of antigenic variation. Agur et al. [3J have studied a model of antigenic 
variation of African trypanosomes which suggests that sequential appearance of different antigenic 
variants can be explained by fitness differences between single- and double-expressors - antigenic 
variants that express one or two VSGs. However, this idea is not supported by the experimental 
evidence arising from normal in vivo growth and reduced immunogenicity of artificially created 
double expressors [17]. Frank (23J El] has suggested a model that highlights the importance of 
cross-reactivity between antigenic variants in facilitating optimal switching pattern that provides 
sequential dominance and extended infection. Antia et al. [1] have considered variant-transcending 
immunity as a basis for competition between variants, which can promote oscillatory behaviour, but 
this failed to induce sequential expression. Many other mathematical models of antigenic variation 
have been proposed and studied in the literature, but the discussion of their individuals merits and 
limitations is beyond the scope of this work. 

In this paper we concentrate a model proposed by Recker et al. [52] (to be referred to as 
Recker model), which postulates that in addition to a highly variant-specific immune response, 
the dynamics of each variant is also affected by cross-reactive immune responses against a set of 
epitopes not unique to this variant. This assumption implies that each antigenic variant experiences 
two types of immune responses: a long-lasting immune response against epitopes unique to it, and 
a transient immune response against epitopes that it shares with other variants. The main impact 
of this model lies in its ability to explain a sequential appearance of antigenic variants purely on the 
basis of cross-reactive inhibitory immune responses between variants sharing some of their epitopes, 
without the need to resort to variable switch rates or growth rates (see Gupta [33J for a discussion 



2 



of several clinical studies in Ghana, Kenya and India, which support this theory). 

From mathematical perspective, certain understanding has been achieved of various types of 
dynamics that can be obtained in the Recker model. In the case when long-lasting immune re- 
sponses do not decay, numerical simulations in the original paper |52| showed that eventually all 
antigenic variants will be cleared by the immune system, with specific immune responses reaching 
protective levels preventing each of the variants from showing up again. Blyuss and Gupta [10J have 
demonstrated that the sequential appearance of parasitemia peaks during such immune clearance 
can be explained by the existence of a hypersurface of equilibria in the phase space of the sys- 
tem, with individual trajectories approaching this hypersurface and then being pushed away along 
stable/unstable manifolds of the saddle- cent res lying on the hypersurface. They also numerically 
analysed robustness of synchronization between individual variants. Under assumption of perfect 
synchrony, when all variants are identical to each other, Recker and Gupta [33] have analysed 
peak dynamics and threshold for chronicity, while Mitchell and Carr |45j have investigated the 
additional effect of time delay in the development of immune response. De Leenheer and Pilyugin 
|18| have replaced linear growth of antigenic variants in the original model by the logistic growth, 
and have studied the effects of various types of cross-reactivity on the dynamics, ranging from no 
cross-reactivity to partial and complete cross-immunty. Mitchell and Carr [JH] have studied the 
appearance of synchronous and asynchronous oscillations in the case of global coupling between 
variants (referred to as "perfect cross immunity" in |18j). 

So far, mathematical analyses of the Recker model have concentrated primarily on identifying 
and studying different types of behaviour in the model. Whilst this has given a certain headway in 
the understanding of possible dynamics, symmetric properties of the model have remained largely 
unstudied, and yet they can provide important insights into the dynamics of the model allowing one 
to distinguish the interactions between variants and the immune system from the effects of topology 
of coupling between antigenic variants. The importance of this topology has been highlighted in 
recent works \12\ [TBI 1ST)] . In this paper we study Recker model from the perspective of symmetric 
dynamical systems. This allows us to perform a systematic analysis of steady states and their 
stability, as well as to classify various periodic behaviours in terms of their symmetries. The outline 
of the paper is as follows. In the next section we formalize the Recker model and discuss some of 
its properties. Section 3 reviews some concepts and techniques from equivariant bifurcation theory. 
In Section 4 we analyse steady states and their stability with account for symmetry properties of 
the system. Section 5 is devoted to symmetry-based classification of different dynamical regimes. 
The paper concludes in Section 6 with discussion of results. 

2 Mathematical model 

Following Recker et al. [52j, we assume that within a human host, the parasite population of P. 
falciparum consists of iV distinct antigenic variants, with each antigenic variant i, 1 < i < N, 
containing a single unique major epitope that elicits a long-lived (specific) immune response, and 
also several minor epitopes that are not unique to the variant. Assuming that all variants have the 
same net growth rate <ft, their temporal dynamics is described by the equation 



where a and a' denote the rates of variant destruction by the long-lasting immune response z% 
(specific to variant i) and by the transient immune response Wi, respectively. The dynamics of the 
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Figure 1: Interaction of malaria variants in the case of two minor epitopes with two variants in 
each epitope. 



variant-specific immune response can be written in its simplest form as 



dzj 
dt 



PVi ~ VZi, 



(2) 



with P being the proliferation rate and \i being the decay rate of the specific immune response. 
Finally, the transient (cross-reactive) immune response can be described by the minor modification 
of the above equation ([2]): 

= P'^Zvj -fJ-'wi, (3) 

where the sum is taken over all variants j sharing the epitopes with the variant i. We shall use 
the terms long-lasting and specific immune response interchangeably, likewise for transient and 
cross-reactive. 

The above system can be formalized with the help of adjacency or connectivity matrix A 
whose entries Aij are equal to one if the variants i and j share some of their minor epitopes and 
equal to zero otherwise |10t I18j . Obviously, the matrix A is always a symmetric matrix. Prior 
to constructing this matrix it is important to introduce an ordering of the variants according to 
their epitopes. Whilst this choice is pretty arbitrary, it has to be fixed before the analysis can be 
done. Consider a system of antigenic variants with just two minor epitopes and two variants in 
each epitope. In this case, the total number of variants is four, and we enumerate them as follows 



1 

2 
3 
4 



11 
12 
22 
21 



(4) 



After the ordering of variants has been fixed, it is straightforward to construct the connectivity 
matrix A of variant interactions. For the particular system of variants Q illustrated in Fig. [TJ this 
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matrix has the form 
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(5) 



With the help of lexicographic ordering, one can systematically construct matrix A for an 
arbitrary number of minor epitopes [ID]. For the rest of the paper we will concentrate on the case 
of two minor epitopes, but the results can be generalized to larger systems of antigenic variants. 
Using the connectivity matrix one can rewrite the system ([!])- ^ in a vector form 




F{y,z, w) 




a w 



(6) 



where y = (yi, y 2 , Un), z = (zi, z 2 , z N ), w = (wi, w 2 , tojv), In denotes a vector of the 
length iV with all components equal to one, and in the right-hand side of the first equation multi- 
plication is taken to be entry-wise, so that the output is a vector again. The above system has to 
be augmented by appropriate initial conditions, which are taken to be 

y(0) > 0,z(0) > 0,w > 0. 

As it has been shown in |1U] . with these initial conditions the system ([6]) is well-posed, so its 
solutions remain non-negative for all time. 

We will assume that cross-reactive immune responses develop at a slower rate than specific 
immune responses, have a shorter life time, and are less efficient in destroying the infection. This 
implies the following biologically realistic relations between the system parameters 



a < a, (j, < fjf } < (3. 



(7) 



3 Elements of equivariant bifurcation theory 

Before proceeding with the analysis of symmetry effects on the dynamics of systems of antigenic 
variants, we recall some concepts and results from equivariant bifurcation theory \30\ 131] . Let 
r C GL(R N ) be a compact Lie group acting on Mr. We say that a system of ODEs 

x = F(x), x£R N , F:R N ^R N , (8) 

is equivariant with respect to a symmetry group T if the vector field F commutes with the action 
of r, i.e. if it satisfies an equivariance condition 

F(jx) = rF(x) for all x G R N , 7 G T. 

The main examples of the symmetry groups we are interested in are T = S n , Z n and D n (of order 
n\, n, and 2n, respectively). Here, S n denotes the symmetric group of all permutations in a network 
with an all-to-all coupling, i.e. all permutations on n symbols. The cyclic group Z n describes the 
symmetry of a unidirectional ring (rotations only), while the dihedral group D n corresponds to a 
bidirectional ring of n coupled units (rotations and reflections in the plane that preserve a regular 
n-gon) . 

One can define the group orbit of a point x G R N under the action of T is defined as 

Tx = {jx G R N : 7 G T} . 
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Note that the equivariance of the system ([8]) implies that its equilibria come in group orbits, since 
if F(x) = for some x G M N , then F(jx) = 7O = for all 7 G T. The isotropy subgroup E(ar) C T 
of a point x G is defined as the subgroup of T which fixes the point x, i.e. 

T,(x) = {7 G T : 7X = x} . 

Associated with each isotropy subgroup S C T is the fixed-point subspace denoted by Fix(S), which 
is the set of points x G l w invariant under the action of S: 

Fix(E) = {x G M. N : ax = x for all a G E} . 

An important property of fixed-point subspaces is that they are flow- invariant, since if x G Fix(E), 
then o~F{x) = F(ax) = F(x), and therefore, F(x) G Fix(E). The equivariant branching lemma 
states that provided certain (generic) conditions are satisfied by the bifurcation, there exists 
a branch of equilibrium solutions with symmetry £ for each isotropy subgroup E C T with 
dim(Fix(£)) = 1. Such isotropy subgroups with one-dimensional fixed-point subspaces are called 
axial. Similarly, one can define C- axial subgroups as those subgroups E C T x S 1 , for which £ 
is an isotropy subgroup of the action of T x S 1 on the centre subspace of the equilibrium, and 
dim(Fix(£)) = 2. 

Next, we recall that a subspace V of M. N is called T -invariant if 

-/V C V, for all 7 G T. 

Two T-invariant subspaces V and W are called Y -isomorphic if there exists a linear isomorphism 
T : V -> W, such that 

T(^v) = 7(Tu), for all v G V, 7 G T. 

A T-invariant subspace V is T-irreducible if the only T-invariant subspaces of V are {0} and V. 
T-irreducible subspaces can be used to perform an efficient decomposition of the phase space 
that would allow block-diagonalization of the linearization matrix, thus simplifying the analysis 
of stability of steady states. The isotypic decomposition proceeds by decomposing M. into T- 
irreducible subspaces Vj so that Mr = Vq © V% ® ■ ■ ■ ® V m . The isotypic components are then 
formed by combining the irreducible subspaces that are r-isomorphic. The isotypic decomposition 
is R N = W ffi W\ ® ■ ■ ■ ® Wk, k <m, where the Wj are uniquely defined |5T] . 

Now, if x(t) is a T-periodic solution of a T-equivariant system Q, then jx(t) is another T- 
periodic solution of Q for any 7 G Y. Uniqueness of solutions implies that either x(t) and jx(t) 
are identical, or there exists a phase shift 9 G S 1 = M/Z = [0, T), such that 

7x(t) = x(t - 6). 

The pair (7, 6) is called a spatio-temporal symmetry of the solution and the collection of all 
spatio-temporal symmetries of x(t) forms a subgroup AcTxS 1 . One can identify A with a pair 
of subgroups, H and K, such that K C H C T. Define 

H = {7 G r : 7{x(t)} = {x(t)}} spatio-temporal symmetries, 
K = {7 G r : jx(t) = x(t) \/t} spatial symmetries. 

Here, K consists of the symmetries that fix x(t) at each point in time, while H consists of the 
symmetries that fix the entire trajectory. Let 

L K = U 7e/AX Fix( 7 ), (9) 
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and let Nr(K) denote the normalizer of K in T: 

N r (K) = {7 G r : = K}. 

The following theorem gives the necessary and sufficient conditions for H and K to characterize 
spatio-temporal symmetries of a periodic orbit. 

Theorem (H/K Theorem |144l30j). Let T be a finite group acting on M. N . There is a periodic solu- 
tion to some T-equivariant systems of ODEs on M> N with spatial symmetries K and spatio-temporal 
symmetries H if and only if 

(a) H/K is cyclic. 

(b) K is an isotropy subgroup. 

(c) dim Fix(if) > 2. //dim Fix(K) = 2, then either H = K or H = N V (K). 

(d) H fixes a connected component ofFix(K)\Lx. 

Moreover, when these conditions hold, there exists a smooth T-equivariant vector field with an 
asymptotically stable limit cycle with the desired symmetries. 



The H/K theorem was originally derived in the context of equivariant dynamical systems by 
Buono and Golubitsky [H] , and it has subsequently been used to classify various types of periodic 
behaviours in systems with symmetry that arise in a number of contexts, from speciation [59] to 
animal gaits and vestibular system of vertebrates |28| . 

Now we can proceed with the analysis of symmetry properties of the system ^ as represented 
by its adjacency matrix A. In the case of two minor epitopes with m variants in the first epitope 
and n variants in the second, the system ^ is equivariant with respect to the following symmetry 
group [TO] 

p _ J ^ m X 171 7^ n ' QqA 

\ S m x S m x Z 2 , m = n. 

This construction can be generalized in a straightforward way to a larger number of minor epitopes. 
It is noteworthy that while within each stratum we have an all-to-all coupling, the full system does 
not possess this symmetry. System ^ provides an interesting example of a linear coupling, which 
does not reduce to known symmetric configurations, such as diffusive, star or all-to-all [50J. A 
really important aspect is that two antigenic systems with the same total number of variants N 
may have different symmetry properties as described by the group T depending on m and n, such 
that iV = ran. The simplest example of different kinds of splitting is given by iV = 12, which can 
be represented asA r = 2x6orasA r = 3x4. 



4 Symmetry analysis of steady states 

In the particular case of non-decaying specific immune response (// = 0), equilibria of the system 
Q are not isolated but rather form an A^-dimensional hypersurface -£?o = {(y> z > w ) £ ^ N '■ y 
w = N } in the phase space [10J. This hypersurface consists of saddles and stable nodes, and in 
addition to the original symmetry of the system it possesses an additional translational symmetry 
along the z axes. The existence of this hypersurface of equilibria in the phase space leads to a 
particular behaviour of phase trajectories, which mimics the occurrence of sequential parasitimea 
peaks in the immune dynamics of malaria [1U\ 152] . 
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Figure 2: (a) Symmetries of the square, (b) Lattice of subgroups of D4 symmetry group. 

When fi > 0, the structure of the phase space of the system ^ and its steady states is drastically 
different. Now, the only symmetry present is the original symmetry T, and the hypersurface of 
equilibria Hq disintegrates into just two distinct points: the origin O, which is always a saddle, and 
the fully symmetric equilibrium 

E = (Y1 N ,Z1 N ,W1 N ), where 

Y = ^ z = w = <P^n c P' 

a/3// + a'n c f3' /i' afi[i' + a'n c /3' fi' a(3fi' + a'n c f3' \x 

Here n c is the total number of connections for each antigenic variant. It has been previously found 
by means of numerically computing eigenvalues of the Jacobian that the fully symmetric equilibrium 
E may undergo Hopf bifurcation as the parameters are varied [10] . It is worth noting that if one 
assumes all variants to be exactly the same, the original system ^ collapses to a system with just 3 
dimensions, and in this case it is possible to show analytically that the fully symmetric equilibrium 
is always stable [S]. This implies that Hopf bifurcation of the fully symmetric equilibrium takes 
place outside the hyperplane of complete synchrony. 

In order to find analytically the boundary of Hopf bifurcation in terms of system parameters, 
as well as to understand the structure of the bifurcating solution, we concentrate on a specific 
connectivity matrix A that corresponds to a particular case of two epitopes with two variants in 
each epitope, as described by the system Q. Despite its simplicity, this case is still important for 
the following two reasons. Firstly, this is the simplest non-trivial combination of antigenic variants 
which can produce interesting dynamics of interactions with the immune system. Secondly, it can 
be used as a paradigm for many two-locus two-allele models [3U |55] . 

In the case of two epitopes with two variants in each epitope, we have m = n = 2 and N = 4, 
and the system (J™ is equivariant under the action of a dihedral group D4, which is an 8-dimensional 



8 



symmetry group of a square. This group can be written as D4 = {1, (, Q 2 , ( 3 , k, k(, k( 2 , kC 3 }, an( A 
it is generated by a four-cycle £ corresponding to counterclockwise rotation by tt/2, and a flip k, 
whose line of reflection connects diagonally opposite corners of the square, see Fig. [2] The group D4 
has eight different subgroups (up to conjugacy): 1, Z4, and D4, as well as = {1, k} generated 
by a reflection across a diagonal, Df = {1,k(} generated by a reflection across a vertical, D?? = 
{1,C 2 ,^,^C 2 } generated by reflections across both diagonals, and D| = {1, C 2 , K C, K C 3 } generated 
by the horizontal and vertical reflections. Finally, the group Z2 is generated by rotation by tt. The 
lattice of these subgroups is shown in Fig. [2|b). D4 has two other subgroups Z2(k( 2 ) = {1, ftC 2 } 
and Z2(«C 3 ) = {1> K C 3 }, which will be omitted as they are conjugate to and Df, respectively. 
There is a certain variation in the literature regarding the notation for subgroups of D4, and we 
are using the convention adopted in Golubitsky and Stewart [30J, c.f. (BJHUEI]. 

The group D4 has four one-dimensional irreducible representations |20t [29] . Equivariant Hopf 
Theorem |30[ I3"T] states that under certain genericity hypotheses, there exists a branch of small- 
amplitude periodic solutions corresponding to each C-axial subgroup T x S 1 acting on the centre 
subspace of the equilibrium. To find out what type of periodic solution the fully symmetric steady 
state will actually bifurcate to, we can use the subspaces associated with the above-mentioned 
one-dimensional irreducible representations to perform an isotypic decomposition of the full phase 
as follows 0162]: 



space 



>12 



... - ,i 3 (i, 1, 1, 1) e m 3 (i, -1, 1, -1) e m 3 (i, 0, -1, 0) e m 3 (o, i, o, -1) 

Jacobian of the linearization of system ^ near any steady state 
is given by 



(12) 



J(S) 
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^eady state E, this Jacobian takes the block form 
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(13) 



where O4 and I4 are 4x4 zero and unit matrices, and A is the connectivity matrix Rather 
than compute stability eigenvalues from this 12 x 12 matrix, we use the isotypic decomposition (12) 
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to rewrite this Jacobian in the block-diagonal form \29\ [62 

J(E) 



( C + 2D 
3 



S 



3 3 3 \ 
C -2D 3 3 
3 C 3 



(14) 



where 
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(15) 



Here, matrix C is associated with self-coupling, and D is associated with nearest-neighbour cou- 
pling. Stability changes in the C + 2D, C — 2D and C matrices describe a bifurcation of the fully 
symmetric steady state E in the even, odd, and V4 subspaces, respectively [62J. Prior to performing 
stability analysis, we recall the Routh-Hurwitz criterion, which states that all roots of the equation 

A 3 + a 1 X 2 + a 2 A + a 3 = 0, 

are contained in the left complex half-plane (i.e. have negative real part), provided the following 
conditions hold [38] 

04 > 0, i = 1, 2, 3, 

(16) 

aia 2 > a^. 

If the last condition is violated, then the above cubic equation has a pair of purely imaginary 
complex conjugate eigenvalues when 



a,i > 0, i 
aia 2 = a 3 , 



1,2,3, 



(17) 



as discussed in Farkas and Simon |2~T1. 



Proposition 1. The fully symmetric steady state E is stable for a' < a' H , unstable for a' > a' H , 
and undergoes a Hopf bifurcation in the odd subspace at a' = a' H , where 



(18) 



Proof. Stability of the fully symmetric steady state E changes when one of the eigenvalues of 



the Jacobian (14) goes through zero along the real axis or a pair of complex conjugate eigenvalues 
crosses the imaginary axis. Due to the block-diagonal form of the Jacobian it suffices to consider 
separately possible bifurcations in the matrices C, C ± 2D. 



For the matrix C given in (15), the characteristic equation takes the form 

A 3 + aiA 2 + a 2 A + a 3 = 0, 

with 

a% = fj, + //, a 2 = a/3 + a'f3' + /x// , a 3 = a(3fi' + a' [3'/j,. 
Clearly, in this case 01,2,3 > 0, and also 

aia 2 — a 3 = af3/i + a'/3'/i' + + /u') > 0, 
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which according to the Routh-Hurwitz conditions ( 16 ) implies that the eigenvalues of the matrix 
C are contained in the left complex half-plane for any values of system parameters. This means 
that the steady state E is stable in the V4 subspace. 

Matrix C + 2D is equivalent to C upon replacing ft' with 3/3', which allows one to conclude that 
the steady state E is also stable in the even subspace for any values of the system parameters. 

Finally, for the matrix C — 2D, the coefficients of the characteristic equation are 

a% = \x + //, d2 = a(3 — a' /3' + ////, 03 = a/3fj/ — ol $ \x. 



Due to biological restrictions on parameters (I7J), it follows that 01,2,3 > 0. Computing a\a,2 — aj, 
gives 

0102 — 03 = a/3/jL + fJ-fJ, ([J, + //) — a'P'fjf. 

When a' < a' H , where 

«H = W , » (19) 



we have 0102 — 03 > 0, which according to (16) implies that the steady state E is stable in the odd 



subspace. When a' = a' H , one has a\a2 = 03, which coincides with the condition (17). Hence, we 
conclude that at a 1 = a' H , the steady state E undergoes a Hopf bifurcation in the odd subspace. 
For a' > a' H , the steady state E is unstable in the odd subspace. ■ 



The implication of the fact that the Hopf bifurcation can only occur in the odd subspace of 
the phase space |62j is that in the system ^ the fully symmetric state E can only bifurcate to an 
odd periodic orbit, for which variants 1 and 3 are synchronized and half a period out-of-phase with 
variants 2 and 4, i.e. each variant is ir out of phase with its nearest neighbours. 

Besides the origin O and the fully symmetric equilibrium E, the system ^ possesses 14 more 
steady states characterized by a different number of non-zero variants y (for a general system with 
N antigenic variants there would be 2 N — 2 of such steady states). To comprehensively analyse 
these steady states and their stability, let us first introduce auxiliary quantities 

Y _ Y = 

1 P'o/ii + aPn'' 2 2P'a'n + aPn n 

and 

Z X ,2 = 0Y lj2 /n, and W lj2 = pY^/n'. 

There are four distinct steady states with a single non-zero variant yi, which all have the isotropy 
subgroup or its conjugate. A representative steady state of this kind is 

E l = (^,0,0,0,^,0,0,0,^,^1,0,^). (20) 

Other steady states E 2 , E% and E4 are related to E\ through elements of a subgroup of rotations Z4. 



Proposition 2. All steady states E\, E 2 , E%, E4 with one non-zero variant are unstable. 



Proof. As it has already been explained, the steady states -©1,2,3,4 all lie on the same group orbit. 
In the light of equivariance of the system, this implies that all these states have the same stability 
type, and therefore it is sufficient to consider just one of them, for example, E\. Substituting the 
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values of variables in E\ into the Jacobian (13) gives the characteristic equation for eigenvalues 
that can be factorized as follows 



(A-0) • (X + fif(X + n') 3 (X + a'Wi-(/)) 2 x 

[A 3 + A 2 Gu + fj!) + A (nfj,' + + a/3')) + l^i + a'/ 3 '/")] = 0. 

It follows from this characteristic equation that one of the eigenvalues is A = 4> > for any values 
of system parameters, which implies that the steady state E\ is unstable, and the same conclusion 
holds for E 2 , E3 and £4. ■ 

Before moving to the case of two non-zero variants, it is worth noting that the symmetry group 
D4 is an example of a more general dihedral group D n of order 2n, for which there are two distinct 
options in terms of conjugacies of reflections. If n is odd, all reflections are conjugate to each other 
by a rotation. However, when n is even (as is the case for D4), reflections split into two different 
conjugacy classes: one class contains reflections along axes connecting vertices, and another class 
contains reflections connecting the sides. These two conjugacy classes are related by an outer 
automorphism, which can be represented as a rotation through n/n, which is a half of the minimal 
rotation in the dihedral group D n [31] . For the D4 group, these two conjugacy classes are reflections 
along the diagonals, and reflections along horizontal/vertical axes. 

Now we consider the case of two non-zero variants, for which there are exactly six different steady 
states. The steady states with non-zero variants being nearest neighbours on the diagramme Q, 
i.e. (1,2), (2,3), (3,4) and (1,4), form one cluster: 

E 12 = (Y 2 , Y 2 , 0, 0, Z 2 , Z 2 , 0, 0, 2W 2 , 2W 2 , W 2 , W 2 ), 
E 23 = {0,Y 2 ,Y 2 ,0,0,Z 2 ,Z 2 ,0,W 2 ,2W 2 ,2W 2 ,W 2 ), 
E M = {0,0,Y 2 ,Y 2 ,0,0,Z 2 ,Z 2 ,W 2 ,W 2 ,2W 2 ,2W 2 ), 
E u = (Y 2 , 0, 0, Y 2 , Z 2 , 0, 0, Z 2 , 2W 2 , W 2 , W 2 , 2W 2 ), 

while the steady states with non-zero variants lying across each other on the diagonals, i.e. (1,3) 
and (2,3) are in another cluster 

£13 = (Y 1 ,0,Y U 0,Z 1 ,0,Z 1 ,0,W 1 ,2W 1 ,W U 2W 1 ), 
E 24 = (0,Y 1 ,0,YiAZiAZ 1 ,2W 1 ,W 1 ,2W 1 ,W 1 ), 

The difference between these two clusters of steady states is in the above-mentioned conjugacy 
classes of their isotropy subgroups: the isotropy subgroup of the first cluster belongs to a con- 
jugacy class of reflections along the horizontal/vertical axes, with a centralizer given by D|, and 
the isotropy subgroup of the second cluster belongs to a conjugacy class of reflections along the 
diagonals, with a centralizer given by D^. 

Proposition 3. All steady states E\ 2 , E 2 %, E34, E14, and also E13 and E 24 , with two non-zero 
variants are unstable. 

Proof. Using the same approach as in Proposition 2, due to equivariance of the system and the 
fact that within each cluster all the steady states lie on the same group orbit, it follows that for the 
analysis of stability of these steady states it is sufficient to consider one representative from each 
cluster, for instance, E\ 2 and E13. 
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Substituting the values of variables in £12 into the Jacobian ( 13 ), one can find the characteristic 
equation for eigenvalues in the form 

(A + a'W 2 - 4>f (A + y'f{\ + ^) 2 (A 2 + fj,\ + a(3Y 2 ) x 

[A 3 + A 2 (/i + y!) + A {n/jf + Y 2 (af3 + 2a' (3')) + Y 2 {apfi' + 2a! P' y)} = 0. 
It follows that this characteristic equation has among its roots an eigenvalue 

A = <p - a W 2 = — — — - — — > 0. 

a/V + 2a' /3'fi 

Since this eigenvalue is positive for any values of parameters, we conclude that the steady state Ei 2 
(and also the steady states E 2 ^, E34, En) is unstable. 

In a similar way, the characteristic equation for the steady state E13 can be written as follows 

(A + 2a'W l - <j>) 2 (A + y) 2 {\ + y') 2 x 

[A 3 + A 2 (/i + //) + A (//// + Yx{ap + a'P')) + Yi(aPfjf + a'p'fi)] 2 = 0. 
Hence, one of the eigenvalues is 

>(a/V - a'P'11) 



X = <j) - 2dWx 



aPjx' + a'P'n 



which is always positive for biologically realistic restrictions on parameters ([7]). Therefore, we con- 
clude that E13 are £"24 are always unstable. ■ 

For three non-zero variants, we again have four different steady states, which have an isotropy 
subgroup or its conjugate. Introducing relative efficacies of the specific and cross-reactive 
immune responses E z and E w [35j 05] : 

aP a'P' 

-Cj?. — , Him - 



one can write the values of state variables for these states in terms of E z and E w as follows: 

E z — E w Ij 2 

yi = — f y > 2/2 = 2/4 = Y, y s = 0, Y 



E z ' yz ^ ' yJ ' E 2 + 2E Z E W -E 2 ' 

Zi = PEz-Ev^ Py Z3 = Q) (21) 

/j, E z % y 

P' ( E z -E w \ & ( E z -E w \ P' 

Wl = ^[2+ y, w 2 = w A = ^(i+ z y, w 3 = 2^y. 

Here, we have chosen the three non-zero variants to be variants 1, 2 and 4. The other three steady 
states can be obtained from this one by rotation of variants. 

Proposition 4. All steady states with three non-zero variants are unstable. 
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0.2 0.4 0.6 0.8 1 

Figure 3: Lines of equilibria in the system ^ for E z = E w . Each vertical line is parameterized by 
s. Parameter values are (f> = 7.5, a = 15, j3 = /?' = 10, // = 0.5 



Proof. Similar to the proofs of Proposition 2 and 3, we employ system equivariance and the fact 
that the steady states with three non-zero variants all lie on the same group orbit to conclude that 



they all have the same stability type. Substituting the values of variables in the steady state (21) 



into the Jacobian (13) yields the characteristic equation 

fx - + 2^^) (A + M )(A + //)P 9 (A) = 0, 

where -Pg(A) is a 9-th degree polynomial in A. It follows that one of the characteristic eigenvalues 

>a 'P'_ 4> («W 2 - a' W) 



2 — -Y 



fj,' a 2 (3 2 (i' 2 - a' 2 f3' 2 fi 2 + 2a(3/ia' f3' >' ' 

which is always positive due to biological restrictions on parameters 0. Hence, the steady state 



(21) and the other three steady states with three non-zero variants are all unstable. 



This completes the analysis of the steady states of system ^ with D4 symmetry for generic 
values of parameters. It is worth noting, however, that in a particular case when E z = E w , the 
determinant of the system of linear equations that describes the steady states is equal to zero, and 
this results in the existence of a line of equilibria 

(f)fj, 4>u 
Vi =V3 = ~ s > 2/2 = 2/4 = s, < s < 



2aj3 ' "* " ' 2a/3 

parameterized by an additional free parameter s. This line is reminiscent of the hypersurface of 
equilibria in the case [i = in that every point on each such line is a steady state of the system ([6]), 
and as one moves along the line of equilibria, stability of the steady states can change, as illustrated 
in Fig. [3} One can observe that as the decay rate of the specific immune responses increases, this 
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Figure 4: Temporal dynamics of the system ([6]). Parameter values are (f> = 7.5, (3 = 10, j3' = 9, 
a = 15, fjL = 0.01, fi' = 0.02. (a) Stable fully symmetric equilibrium (a' = 8). (b) Anti-phase 
periodic solution (a' = 8.5), spatio-temporal symmetry (H, K) = (D4, D|). (c) Anti-phase periodic 
solution (a' = 11.6), spatio-temporal symmetry (H,K) = (D4, Dg). (d) Discrete rotating wave 
(a' = 11.6), spatio-temporal symmetry (H,K) = (Z4, 1). (e) Out-of-phase oscillations (a' = 11.6), 
spatio-temporal symmetry (H,K) = (D^D^). (f) Chaos (a' = 12). 



leads to stabilization of such steady states, and provided this rate is sufficiently high, all steady 
states on lines of equilibria are stable. 

5 Dynamical behaviour of the model 

In the previous section we established that a fully symmetric steady state E can undergo Hopf 
bifurcation, giving rise to a stable anti-phase periodic solution. Now we look at the evolution of 
this solution and its symmetries under changes in system parameters. For convenience, we fix all 
the parameters except for the rate of variant destruction by transient immune response a' , which 
is taken to be a control parameter. The results of numerical simulations are presented in Fig. |4j 
When a' is sufficiently small to satisfy the condition a' < a' H , the fully symmetric steady state is 
stable, as shown in plot (a). As it crosses the boundary E w = E z , the fully symmetric steady state 
loses stability through a Hopf bifurcation giving rise to a periodic solution with the spatio-temporal 
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symmetry H = D4 and a spatial symmetry K = D|, which ensures that the variants 1 and 3 have 
the same behaviour, as do variants 2 and 4, as is illustrated in figure (b). Figure (c) indicates that 
as a' increases, the periodic solution retains its symmetry but changes temporal profile, acquiring 
two peaks during a single period. 

When analysing other types of periodic behaviour in the model, one can note that according 
to the H/K theorem (see Section [3]), periodic states can have spatio-temporal symmetry group 
pairs (H,K) only if H/K is cyclic, and K is an isotropy subgroup [TOj, [30]. In the case of D4 
symmetry group acting on four elements, there are eleven pairs of subgroup H and K satisfying 
these requirements [30]. Now we look at actual solutions of the model and identify these spatial 
and spatio-temporal symmetries. For the same value of a' as in Fig. (c), the system exhibits 
two more solutions with quite different spatial and spatio-temporal symmetries, as demonstrated 
in Figs, (d) and (e). The solution shown in Fig. (d) has the symmetry (H,K) = (Z4, 1) and 
is a discrete travelling wave, also known as a "splay state" [61], "periodic travelling (or rotating) 
wave" [7], or "ponies on a merry-go-round" or POMs [5] in the studies of systems of coupled 
oscillators. In this dynamical regime all variants appear sequentially one after another along the 
diagramme ([!]) with quarter of a period difference between two neighbouring variants. From the 
perspective of equvariant bifurcation theory, this solution is generic since the group Z n is always 
one of the subgroups of the D n group for the ring coupling, or the S n group for an all-to-all 
coupling, and its existence has already been extensively studied [5j [29l [3T] . From the immunological 
point of view, this is an extremely important observation that effectively such solution, which 
immunologically represents sequential appearance of parasitemia peaks corresponding to different 
antigenic variants, owes its existence not to the individual dynamics of antigenic variants, but rather 
to the particular symmetric nature of cross-reactive interactions between them. This immunological 
genericity ensures that the same conclusions hold for a wide variety of immune interactions between 
human host and parasites, which use antigenic variation as a mechanism of immune escape, as 
illustrated, for instance, by malaria, African Trypanosomes, several members of Neisseria family 
(N. meningitidis and N. gonorrhoeae), Borrelia hermsii etc. [33^ 163] . 

Another co-existing solution for the same value of a' is a state shown in Fig. (e). This solutions 
is characterized by the symmetry (H, K) = (Dg, D?) and corresponds to a situation, in which 
variants 1 and 3 are oscillating half a period out-of-phase with each other, and the variants 2 and 
4 coincide and oscillate at twice the frequency of the pair (1,3). As the value of a' is increased 
further, the system demonstrates a state of deterministic chaos, when variants appear in arbitrary 
order and magnitude. These kinds of periodic and chaotic solutions have been previously identified 
as most relevant for the clinical analysis of blood-stage malarial infection [57J. 

6 Discussion 

In this paper we have used techniques of equivariant bifurcation theory to perform a comprehensive 
analysis of steady states and periodic solutions in a model of antigenic variation in malaria. In 
the simplest case of two epitopes with two variants in each epitope, the system is equivariant with 
respect to a D4 symmetry group of the square. Using isotypic decomposition of the phase space 
based on the irreducible representations of this symmetry group has allowed us to find a closed 
form expression for the boundary of Hopf bifurcation of the fully symmetric steady state in terms 
of system parameters, as well as to show that for any parameter values this bifurcation leads to an 
anti-phase periodic solution. All other steady states have been classified in terms of their isotropy 
subgroups, and this has provided insights into their origin and relatedness. With the help of the 
H/K theorem, we have classified various periodic states of the model in terms of their spatial and 
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temporal symmetries. 

One of the important issues that have to be taken into account when applying methods of 
equivariant bifurcation theory to the models of realistic systems is the fact that these systems do 
not always fully preserve the assumed symmetry. In the context of modelling immune interactions 
between distinct antigenic variants, this means that not all variants cross-react in exactly the 
same quantitative manner. Despite this limitation, due to the normal hyperbolicity, which is a 
generic property in such models, main phenomena associated with the symmetric model survive 
under perturbations, including symmetry-breaking perturbations. The discussion of this issue in 
the context of modelling sympatric speciation using a symmetric model can be found in [30J . 

Whilst the analysis in this paper was constrained to the models of within-host dynamics anti- 
genic variation, a similar approach can also be used for the population level mathematical models 
of multi-strain diseases. These models usually have a similar structure in that all strains/sero types 
are assumed to be identical and have a certain degree of cross-protection or cross-enhancement 
based on antigenic distance between them [D [21 [Tg [161 [23 [32l E3 HH [53l [55]. The fact that mul- 
tiple strains are antigenically related introduces symmetry into a model, and this then transpires 
in different types of periodic behaviours observed in these models. To give one example, almost 
all of these models support a solution in the form of discrete travelling wave having a symmetry 
(H,K) = (Z„,l), which represents sequential appearance of different strains [351 HH [55]. I n a 
model studied by Calvez et al. [15] (based on an earlier model of Gupta et al. [35]), the authors 
study a multi-strain model having symmetry of a cube, and they find that the fully symmetric 
steady state can undergo Hopf bifurcation and give rise to an antiphase solution with a tetrahedral 
symmetry, which was not expected a priori. At the same time, if one looks at this system from 
the symmetry perspective, such periodic solution is to be expected as it has a D4 symmetry, which 
corresponds to one of the three maximal isotropy subgroups and is therefore generically expected 
to arise at a Hopf bifurcation [22] [37] . Observations of this kind illustrate that taking into account 
symmetries of the underlying models of multi-strain diseases can help make significant inroads in 
understanding and classifying possible periodic behaviours in such systems. 

The results presented in this paper are quite generic, and the conclusions we obtained are valid 
for a wide range of mathematical models of antigenic variation. In fact, they are applicable to 
the analysis of within-host dynamics of any parasite, which exhibits similar qualitative features 
of immune interactions based on the degree of relatedness between its antigenic variants. The 
significance of this lies in the possibility to classify expected dynamical regimes of behaviour using 
very generic assumption regarding immune interactions, and they will still hold true provided the 
actual system preserves the underlying symmetries. In the model studied in this paper the degree of 
cross-reactivity between antigenic variants does not vary with the number of epitopes they share. It 
is straightforward, however, to introduce antigenic distance between antigenic variants in a manner 
similar to the Hamming distance [H [5"3"1 158] . Such a modification would not alter the topology 
of the network of immune interactions but rather assign different weights to connections between 
different antigenic variants in such a network. Symmetry analysis of the effects of antigenic distance 
on possible dynamics are the subject of further study. 
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